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ABSTRACT 

Numerical simulation of magnetohydrodynamic (MHD) turbulence makes it 
possible to study accretion dynamics in detail. However, special effort is required 
to connect inflow dynamics (dependent largely on angular momentum transport) 
to radiation (dependent largely on thermodynamics and photon diffusion). To 
this end we extend the flux-conservative, general relativistic MHD code HARM 
from axisymmetry to full 3D. The use of an energy conserving algorithm allows 
the energy dissipated in the course of relativistic accretion to be captured as 
heat. The inclusion of a simple optically thin cooling function permits explicit 
control of the simulated disk's geometric thickness as well as a direct calcula- 
tion of both the amplitude and location of the radiative cooling associated with 
the accretion stresses. Fully relativistic ray-tracing is used to compute the lu- 
minosity received by distant observers. For a disk with aspect ratio H/r ~ 0.1 
accreting onto a black hole with spin parameter a/M = 0.9, we find that there 
is significant dissipation beyond that predicted by the classical Novikov-Thorne 
model. However, much of it occurs deep in the potential, where photon capture 
and gravitational redshifting can strongly limit the net photon energy escaping 
to infinity. In addition, with these parameters and this radiation model, sig- 
nificant thermal and magnetic energy remains with the gas and is accreted by 
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the black hole. In our model, the net luminosity reaching infinity is 6% greater 
than the Novikov-Thorne prediction. If the accreted thermal energy were wholly 
radiated, the total luminosity of the accretion flow would be ~ 20% greater than 
the Novikov-Thorne value. 

Subject headings: Black holes - magnetohydrodynamics - instabilities - stars: accretion 



1. Introduction 

For the past thirty-five years, it has been the standard view in the astrophysical com- 
munity that the total amount of energy per unit mass dissipated in the course of accretion 
onto a black hole is exactl y equal to the binding energy of the innermost stable circular orbit 
( iNovikov fc Thorndll973l ); consequently, it depends only on the black hole spin parameter 
a/M. The argument leading to this result depended on a number of assumptions: that 
the flow is time-steady and axisymmetric, that any heat dissipated is promptly radiated, 
and that the r-0 component of stress goes to zero at the ISCO. Although the first several 
assumptions appear t o be relatively innocuous, the last was regarded as questionable almost 
from the be ginning ( Thorn e 1974) and h as been subject to renewed questioning in more 
recent years (jKrolikill999bl ; lGammidll999l ). If accretion flows were fundamentally hydrody- 
namic, the heuristic argument for this boundary condition (that the inertia of matter in the 
plunging region should always be much smaller than the inertia in the stable-orbit portion 
of the disk) would be cogent; the point raised by all its critics is that if magnetic fields play 
an important role, their stress would not necessarily be diminished even in regions of small 
matter density. Resolution of this point is important because continued forces at and inside 
the ISCO would permit continued dissipation, possibly substantially increasing the total. 

Although the significance of magnetic forces was no more than a speculation when 
the zero-stress boundary condition was first critic ized, in recent years it has been recog- 



nized that they are, in fact, essential to accretion (IBalbus fc Hawleylll998l ). Stimulated by 



this recognition, the past decade has seen many numerical simulations of global disk dy- 



ics (MHD) ( 


Hawlev & Kro 


ik 


2001. 


2002 


: lArmitaee et al. 


2001: 


Reynolds & Armitaee 


2001: 


Armitaee & Reynold 


\ 2003 


Machida k Matsumotoll2003l: 


De Villiers et al. 


2003; 


Krolik et al. 


2005; 


Gammie et al. 


2004) 


Initially these simulations assumed Newtonian dynamics in a 



pseudo- Newtonian potential; in the m iddle of this effort, new codes were developed that per- 
mit simulations in general relativity (jPe Villiers &: Hawleyl 120031 ; iGammie et al.ll2003l ). So 
far, while often treating angular momentum transport quite accurately, all of these simula- 
tions have handled thermodynamics and energy transport comparatively crudely: In GRMHD, 
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the code developed by De Villiers and Hawley, only an internal energy equation is solved, in 
which the gas is assumed to behave adiabatically except in shocks; in this code, therefore, 
there are sizable (and uncontrolled) numerical losses of energy whenever magnetic field or 
kinetic energy is lost on the gridscale. By contrast, in HARM, the code developed by Gam- 
mie et al., a total energy equation is solved (so no energy is lost), but there are also no 
radiative losses. The best effort that could be made to estimate actual radiative efficiency 
was therefore through plausible, but ad hoc models, usually keyed to the magnetic stress 



flBeckwith et al.ll2008a[ ). 



In an effort to remedy this situation, we have altered the HARM code in two significant 
ways. First, we have extended it from 2D (axisymmetric) to 3D. This extension has two 
major consequences: we can study nonaxisymmetric fluctuations, and are free from 2D 
artifacts like the "channel solution" ; and we are not limited by the anti-dynamo theorem to 
short duration simulations. Second, we have introduced a toy-model optically thin cooling 
function. By this means, we can track how much radiation might be produced (and where) 
in order to compute the radiative efficiency explicitly. We can also use this cooling function 
to regulate the geometric thickness of the accretion flow. A detailed description of the new 
code (which we call HARM3D) can be found in § 2. 

For our first use of this code, we chose to run a simulation that would illustrate how 
MHD turbulence influences the global energetics of accretion onto a black hole. Its results 
can be compared directly to those of NT: Time-averages of its data can be matched against 
the classical model's steady state. Quantities integrated over constant-radius shells can be 
compared with the corresponding vertically-integrated ones derived assuming axisymmetric 
and "razor-thin" disks. The cooling function can be designed to (almost) reproduce the 
prompt radiation assumption. However, we have no need to impose any guessed boundary 
condition on the stress because in this numerical calculation we are able to use the the real 
physical boundary condition to accretion dynamics: the black hole's event horizon. Thus, 
the ratio of the energy radiated in this simulation to the mass accreted in it provides a 
direct test of how much the zero-stress boundary condition affects the radiative efficiency. In 
addition, of course, we will also be able to examine the interesting effects of non-stationary 
flow, non-axisymmetry, and so on. 

Because we recognize that quantitative results may well depend on a number of pa- 
rameters (magnetic field configuration and disk thickness, most notably) and because our 
radiation model does not fully represent any particular physical situation, we emphasize that 
the numbers we present here are only preliminary samples. When we discuss these results, 
we will explain more specifically the degree to which they are model-dependent. We intend 
to explore more fully in future work both how to model this process more realistically and 
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how external parameters such as magnetic configuration and accretion rate couple with black 
hole spin to control the radiative output of accretion onto black holes. 



2. The Computation: HARM3D and the Parameters of Our Simulation 

Q uite a number of general relativistic MHD simulation codes have been written al 



readv (Komissarov 
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De Villiers & Hawlev 2003 




Duez et alJ2005 


Shibata & Sekieuchi 


2005; 


Anninos et al. 


2005 


Anton et al. 2006: Noble et al. 


2006: Mizuno et al. 


2006; 


Anderson et al. 


2006; 


Tchekhovskov et al. 2007: Frarile et al.ll2007 




Del Zanna et al. 


20071: Cerda-Duran et a 


. 2008 


). Our starting point for the code used in this 


paper was the HARM code 


Gammie et al. 


2003|; 


Noble et al 


2006 


). HARM solves the equations 



mentioned, axisymmetric calculations suffer from two major draw backs: the dominance o f 



"channel solutions", which are ubiquitous in 2D but unstable in 3D (IBalbus fc Hawleyll 19981 ). 
and the fact that neither turbulence nor magnetic field can be sustained indefinitely in 2D. 
To avoid these limitations, we extended the algorithm to three spatial dimensions. HARM's 
conservative formulation means that it does not lose energy to numerical dissipation; rather, 
kinetic and magnetic energies lost at the gridscale are captured as heat. At the same time, a 
conservative formulation permits easy introduction of a formal radiative cooling term. Thus, 
our new code, called HARM3D, is the first global MHD accretion code to treat thermodynamics 
in a controlled fashion^. 



2.1. Basic Equations 



We begin the descriptio n of HARM3D wit h an ex plicit statem ent of the equations governing 
our model. Contrasts with Gammie et al. J2003h (HARM) and be Villiers fc Hawlevi J2003h 
(GRMHD) will be highlighted along the way. We use Greek letters for spacetime indices, and 
Roman letters for spacelike indices. The signature of the metric is the same as the one used in 



1 Technically, our code's pro cedural flow and data structure design developed from an early 3D version 
of the HAM code ( Gammielbood ) that is now publicly available as a shearing box code (jGammiel 1 1999-20081 



Guan fc Gammie 



200- 



HARM found in lGammiej (jl999-2008h . 



H). All other routines were either developed by us, or taken from the public version of 



2 Concurrent with our effort, IShafee et all (|2008l ) have also developed a 3D version of the HARM code. 
They have now used this to explore the dynamical effects of cooling in a disk around a Schwarzschild black 
hole, but did not directly measure the radiation it produced. 
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Misner et al.l (119701 ) (i.e., -+++), and geometrized units are used such that G — c — M — 1. 

The general relativistic MHD (GRMHD) equations of motion include the continuity 
equation, 

V M (pu") = , (1) 
the equations of local energy conservation 

V^ = , (2) 

and Maxwell's equations 

V^ = , (3) 

V„F^ = J* 1 . (4) 

Here, p is the rest-mass density, w M is the 4-velocity of the fluid, F^ u is the Faraday tensor 
times l/\/47r, F v is the dual of this tensor or the Maxwell tensor times 1/a/47t, and J M is 
the 4-currento The total stress-energy tensor is the sum of the fluid part, 

TZ d = pWV + Pgr, (5) 

and the electromagnetic part 

TZi = F^ X F\ - ^F Xk F Xk = \\b\\W + ^\\b\\ 2 g^ - IfW , (6) 



where g^ v is the metric, h = (1 + e + P/p) is the specific enthalpy, P is the pressure, e is the 
specific internal energy density, = F u ^u u is the magnetic field 4-vector, and ||6|| 2 = b^b^ 
is twice the magnetic pressure P„fl 

Equations (JT}{3]) can be expressed in flux conservative form 

9 i U(P) = -9 4 P(P) + S(P) (7) 

where U is a vector of "conserved" variables, F* are the fluxes, and S is a vector of source 
terms. Explicitly, these are 

U (P) = [pu\ T\ + pu\ T f j, B k ] T (8) 



3 We follow Gammie et al. ( 20031 ) in our definition of the electromagnetic field tensor and magnetic field 
variables. 

4 The magnetic 4-vector 6 M defined in this paper is equivalent to that in HARM and GRMHD, even though our 
and HARM's definition is different from GRMHD's by a sign. For this reason, GRMHD's version of our equation (TJ| 
differs by a sign. These sign differences can all be reconciled by noting that their electromagnetic field tensors 
have opposite sign. The resulting equations of motion are independent of these sign conventions. 
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F* (P) = ^ [pu\ T\ + pu\ T v {Vu k - b k v?)} 
S (P) = [0, T\T\ K , T\T X ]K , 0] T 



(9) 
(10) 



where g is the determinant of the metric, T x tlK is the metric's affine connection, and B % = ~F' 



is our ma. 



gnetic field B Note that the source term for the energy equation is non-zero only 
when the metric is time-dependent (as evidenced by its proportionality to r^J. The equations 
of motion are closed by an equation of state, P — (Y — 1) pe, where T is the adiabatic index, 
set to 5/3 in this work. The primitive variables, P = j/?, P, u 1 }, a r e rec overed using an 
optimized version of the the "2D" algorithm described in lNoble et al.l (120061 ) . The primitive 
velocity is the flow's velocity as viewed by a zero angular momentum observer (ZAMO): 



u* = u { + aWg* 



where a = 1/ y/—g tt is the lapse function and W = cm* is the Lorentz factor. 



2.2. Initial Data 



In the initial state of the simulation, the matter is in an axisymmetric hydrostatic 
torus that orbits the black hole with a specific angular momentum profile slightly shal- 
lower than Keplerian and u r = u e = 0. The disk is centered about the equator of the 
black hole's spin and is initially assumed to be isentropic. In curved spacetimes, the an- 
gular frequency — Q = u^/u 1 — is not a simple function of the specific angular momentum — 
/ = —Utft/ut- For example, one can show that when u r = u 9 = in Boyer-Lindquist coordi- 
nates, £1 = {g t< t > — g^l) I (g u — g t( ^l) . In or der to solve the time -inde pendent Euler equations , 
we must therefore specify l(r, 6). Following IChakrabartil (119851 ) and lDe Villiers et al.l (120031 ). 
we do this by assuming that Q ~ X~ q , where A 2 = l/Q. The solution is simplified by setting 
A to its Schwarzschild value A = 



V 9 U V b 1 i which is exact when a = but leads to a 
solution marginally out of equilibrium when a ^ 0; the slight departure from equilibrium 
insignificantly affects the disk's evolution because the magnetic field quickly becomes dy- 
namically important. Ultimately, we arrive at an equation for l(r,6): l/i m = (\/ \ m ) 2 ~ 9 , 
where l in = /(r in ,7r/2) and Ai n = A(r in ,7r/2). 

With the in t ention of closely mimicking the initial conditions of simulation KDP of 
De Villiers et al.l (120031 ). we put the torus pressure maximum at r = 25M and choose an 
angular momentum distribution parameter q = 1.67. The torus inner boundary is r in = 15M, 



5 The "CT field" of GRMHD, B\ is proportional to our magnetic field: B i = y /-AwgB i 
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with l m = 4.576. These parameters yield a disk very similar to that of De Villiers et al., but 
with a slightly larger l in . 

The solution to Euler's equations provides us with h and u^. The rest-mass density 
is then calculated from the equations of state — P = (T — 1) pe and P = Kp v — and h: 
p = [{h — 1) (r — 1) / (.ftT)]^ 1 '^ -1 . We suppose that the gas is non-relativistic, choosing 
T = 5/3 and K = 0.01. Integrating over the volume of the initial gas distribution, we find a 
total rest-mass of 353. This is 20% larger than that in simulation KDP, a shift due to our 
slightly different choice of l- m . Note that the code units of gas mass are completely arbitrary. 

The initial magnetic field lies entirely within the torus and follows contours of constant 
density. The magnitude of the magnetic field is set so that the volume-weighted integrated 
magnetic pressure is 100 times less than the volume-weighted integrated gas pressure. 

The atmosphere surrounding the disk is unmagnetized and static. The atmosphere's 
density and pressure are set to their smallest allowed values, which are chosen so that the 
floor state is in approximate pressure equilibrium: p Rooj . = 7 x 10~ 9 p max r~ 3 / 2 and Pfl oor = 
7 x 10 _11 p max r~ 5 / 2 (r — 1), where p max is the initial maximum value of the rest-mass density 
in the disk. 



2.3. Radiative Cooling 

A magnetized accretion disk is subject to the magneto-rotational instability (MRI), 
which transfers angular momentum outward. This transfer taps into the available free energy 
of differential rotation, creating the magnetic fields and poloidal velocity fluctuations that 
make up the resulting MHD turbulence. This turbulence is dissipative; magnetic and kinetic 
energy is lost numerically at the gridscale. Equation ([7j), however, ensures that in the 
numerical solution all that dissipated energy is converted to heat. If that heat were retained 
by the fluid, the disk would become ever hotter and geometrically thicker. Ultimately the 
thermal energy would either be accreted by the hole or be carried out from the disk by a 
wind. By adding a loss term to the energy equation, we can estimate either the luminosity 
of those systems in which radiation is efficient or the total heat generated in those systems 
in which it is not. We assume that the radiation described by this loss term is optically thin. 
It therefore acts as a passive sink in the local energy conservation equation (J5J): 

V^ v = -T v , (12) 

where T v is the amount of radiated energy-momentum per unit 4-volume in the coordinate 
frame. To describe the radiation, we make the simplest assumption: that the emission is 
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isotropic in the fluid's frame: 

T v = Cu u (13) 

where the "cooling function" C is the rate energy is radiated per unit proper time in the 
fluid frame. 

The NT assumptions include complete prompt radiation of all locally-dissipated heat. 
We cannot exactly replicate that in a simulation, for the gas must retain some thermal energy. 
However, we can arrange for the great majority of the heat to be radiated by constructing 
a cooling function that keeps the temperature of the gas at a small fraction of the virial 
temperature. In so doing, we can also control the disk's aspect ratio H/r, a parameter often 
considered significant in analytic disk models. 

In different contexts, different definitions of the scale-height H are sometimes used. 
For a thin isothermal disk in a Newtonian potential, the density profile is Gaussian, p oc 
exp[— z 2 /(2H G )], with H G = q/Q, for isothermal sound speed q. Another common measure 
of the scale-height is the half-width at half-maximum (HWHM), H HWHM = y/2W!H G . A 
third is the vertical density moment, 

H = J d0d4y/=9Py/9re\0-'Kl1\l J d9d<p^jp. (14) 

When the profile is Gaussian, H = ^/2/ttH g = 0.798# G . 

We prefer the moment definition because it is a direct measure of the characteristic mass- 
weighted disk thickness, it is robust with respect to fluctuations, and it is closely related to 
the characteristic scalelengths of hydrostatic balance. Ideally, it would be computed in the 
fluid-frame, but in the interest of computational economy we define it in the coordinate frame. 
Moreover, when any of these definitions of disk thickness is taken in ratio to the radius, it 
should be recalled that the radial coordinate r is not a proper distance. Unfortunately, there 
is no obvious adequate substitute. 

In any event, given this definition, the temperature that should produce a desired aspect 
ratio H/r in Newtonian gravity is 



T,(r) = \ 



— rfl(r) 
r 



(15) 



In code units, T* = (r - l)e = (2/3)e. 



In our simulation, we evaluate T* in the disk body using the relativistic orbital frequency 
tt(r > r isco ) = 1/ (r 3 / 2 + a/M). In a more completely relativistic treatment, Q(r) would be 
replaced by QkR^ 2 ^), where Qk is the Newtonian Keplerian rotation frequency and R z is 
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the rel ativistic corre ction factor for the vertical gravity ( flAbramowicz et al.lll997l ); notation 



as in (Krolik 1999a ))FI. Inside the innermost stable circular orbit (ISCO), we define f2 as 



the orbital frequency of a particle with the specific energy and angular momentum of the 
circular orbit at the ISCO: 

n(r < *»> = ^j^m; - (16) 

Here is the 4-velocity of the ISCO orbit. 

To ensure that the disk stays near the target temperature, the cooling rate should be 
rapid (i.e., ~ Q), but drop to zero when the temperature falls below T*(r). All of these 
criteria are satisfied by a cooling function with the form 

C = sQpe[Y- 1 + \Y - \\} q , (17) 

where Y = (T — l)e/T*(r) and s is a constant of proportionality. Note that the term in the 
square brackets serves as a switch, so that £ = whenever Y < 1. The exponent q controls 
how rapidly the cooling rate grows when the temperature exceeds the target. We found that 
q = 1/2 cools the plasma efficiently while maintaining a stable evolution, and we set s = 1. 
Only those fluid elements on bound orbits — where (1 + e + P/p) u t > — 1 — are cooled. 

In addition to controlling the vertical thickness of the disk, the cooling function provides 
a self-consistent way of comparing emission from the simulated disk with that expected in a 
standard NT model. When making this comparison, we use the angle-averaged fluid-frame 
luminosity per unit area (of an annulus located at the equator) from our 3D simulation data: 

J ax w |(9= w /2 

where each component of the vector dx^ = e^\dx u represents the extent of a cell's di- 
mension as measured in the fluid element's rest frame, and is the orthonormal tetrad 
that transforms vectors in the Boyer-Lindquist coordinate frame to the local fluid frame (see 



Beckwith et al.l (j2008al ) for explicit expressions for the tetrad) . The vector dx v is the Boyer- 
Lindquist coordinate frame version of the Kerr-Schild vector dx^s = [0, Ar, A8, A<p] (r, 6, (ft), 
where Ar, A8, A(f> are the radial, poloidal and azimuthal extents of our simulation's finite 
volume cell located at (r, 9,cf)). 

We also wish to calculate the radiated luminosity measured by a distant observer in order 
to include the effect of photon losses into the black hole. This is done by ray-tracing through 



6 The expression given in these references contains a typo: should be E^. In addition, we define R z 
inside the ISCO by setting and u t to the values they have at the ISCO. 
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the spacetime and integrating the radiative transfer equation along geodesies. Redshift 
factors from differences in inertial reference frames are automatically taken into account and 
include such effects as gravitational redshift and relativistic beaming. As with the cooling 
function, we assume that the fluid is optically thin and — consequently — ignore scattering 
and absorption. 

To briefly summarize our method, we trace a large number of rays from observers at 
infinity at 8 polar angles and integrate the transfer equation along each ray. Since we aim 
only at estimating the bolometric luminosity of the disk, and not at computing its spectrum, 
we assume that all radiated energy is emitted at a single frequency equal to the Doppler 
shifted frequency of observation. From the transfer solution along each ray, we construct 
images of the disk as it would be seen by each of those observers, and then sum the radiation 
they receive. In order to compute the radiation reaching infinity for the NT model (whose 
photons are also subject to possible capture by the black hole and Doppler shifting), we place 
an emissivity designed to match the NT surface brightness in the two planes of cells nearest 
the equatorial plane. Assuming that the four-velocities of those cells are exactly those of 
circular orbits at those radii, we then compute the luminosity at infinity in this model by 
the same ray-tracing technique as employed on our simulation data. Additional details are 
given in Appendix |A] 



2.4. Coordinates, grid, and boundary conditions 

Equation ((7j) is solved using finite volume techniques on a uniform grid in the so -called 



"Modified Kerr-Schild" (MKS) coordinate system described in iGammie et al.l (120031 ). It is 
based on the Kerr-Schild (KS) coordinate system that eliminates the coordinate singularity 
at the horizon. The modification allows us to adjust the radial and angular discretization 
through a continuous coordinate transformation. We set the MKS parameter huKS = 0.3 



( IGammie et al.ll2003l ). which makes the poloidal cell scale at the axis about 5.7 times larger 
than that at the equator and allows us to resolve greater detail in the accretion disk than 
would be possible with the same number of equally spaced grid cells. 

The simulation reported here used 192 x 192 x 64 cells in the radial, poloidal, and 
azimuthal directions respectively, with r G [1.28M, 120M], 6 G [0.05tt, 0.95tt], <p G [0,tt/2]. 
The radial extent is as large as the one used in KDP except our coordinates penetrate the 
horizon by five cells. We tested whether our polar angle discretization adequately resolved 
the fastest-growing mode of the magneto-rotational instability by calculating — in the local 
fluid frame — the fastest-growing mode's wavelength and the local poloidal size of a cell. 
Averaged over azimuth and time (over t = [7000M, 15000M]), the fastest-growing mode 
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was resolved by at least seven cells at all radii. The absolute minimum number of cells per 
wavelength for all time and radii is never smaller than four. For the time discretization, we 
have found a Courant factor of 0.8 is adequate when used with the existing step size control 
method in HARM. 

In the original MKS coordinate system, cells are placed all the way to the axis. We 
have introduced a new reflecting boundary condition that allows us to excise the coordinate 
singularity there. With the boundary placed at an angle of 0.057T from the axis (as in KDP), 
the excision enlarges the time step we can take, speeding up the evolution by about a factor of 
four relative to simulations using grids without the cut-out. The radial boundary conditions 
are the same as in the released version of HARM, and we use periodic boundary conditions for 
the azimuthal boundaries. 



2.5. Algorithmic details 



The equations of motion are integrated using alm ost the same high-resolution shock- 
capturing methods as described in lGammie et al.l (120031 ) . We still use HARM's Lax-Friedrichs 
numerical flux formula, as it is more diffusive than the HLL formula and seems to be stabler 
for our purposes. How ever, the piecewise-linear re construction is replaced with a parabolic 
interpolation method (IColella &: Woodward! 119841 ) as our means of reconstructing values 
at cell faces. As in HARM, we use an MC (monotonized central-differenced) slope limiter. 
Parabolic reconstruction i mproves stability in low density regions where ||6|| 2 /p >> 1, such 
as are found in the funnel (iMcKinneyl 120061 ) . 



We also use parabolic interpolation in the "Flux-CT" scheme of iTothl (120001 ) that pre- 
serves the divergence constraint. Originally, the electromotive forces (EMFs) at the cell faces 
were calculated by a second-order accurate two-point averaging procedure. This method 
failed to dissipate a cell-scale sawtooth instability seen in the magnetic field along the in- 
tersection between the inner radial and poloidal boundaries. Parabolic interpolation of the 
EMFs, however, is successful at quelling the instability. 

Even with the improvements described so far, stably evolving plasma whose total energy 
is dominated by magnetic and kinetic energies is difficult. In a conservative code like HARM3D, 
the critical step is deriving good primitive variables, P, from the conserved quantites. For 
instance, the magnetic field is typically a few orders of magnitude larger than the pressure 
in the funnel. The pressure is recovered from inverting the equation for the total energy, T* t , 
which involves subtracting T^ Mi and the fluid's inertia term from T l t . In the funnel, this 
operation is essentially a subtraction of two large numbers whose result will likely be the size 
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of either term's truncation error. The subtraction can result in either positive or negative 
pressures. This is known as the "positive-pressure problem" in hydrodyn amics and MHD 
and has been studied extensively (IRyu et al.lll993l ; iBalsara fc Spicerlll999l ). We have found 
that even when positive pressures are recovered, numerical errors may result in pressure 
fluctuations that differ by orders of magnitude between adjacent cells. These fluctuations 
create pressure gradients that can accelerate matter and add energy to the system artificially. 

In order to treat the positive pressure problem and correct for other unphysical states 
that may arise (e.g. W < 0), we have completely redesigned HARM's recovery procedure. The 
most significant change is the inclusion of the conservation of entropy equation 



V M (Sit") = 



where 



S = 



P 



P 



r-i 



(19) 



(20) 



Following a method similar to that of IBalsara &: Spicerl (119991 ). we integrate equation ([191 
in parallel with (J7|). Whenever the standard primitive variable method fails to converge, u l 
is unphysical, or pe < 10~ 2 ||5|| 2 , we use a new inversion method which is identical to the 
standard one except the total energy equation is replaced by equation ([20]) . Even though 
this inversion method is guaranteed to yield a positive pressure, it can either fail to converge 
to a solution or yield an unphysical u l . If either happens, we interpolate P using data from 
neighboring cells for which we have successfully calculated P. Finally, we impose a floor on 
the pressure and density and ensure W < 50 by renormalizing u l . 



We note that using equation ([19]) leads to a method that no longer conserves total 
energy to round-off error, but the impact of these departures from strict conservation is 
limited. The entropy equation is substituted for the energy equation only where the fluid is 
very strongly magnetically-dominated, and only when no energy- conserving method yields a 
physical solution. In the simulation reported here, the net injection or loss of mass, energy 
and angular momentum is only ~ 0.001 — 0.007 times the flux of these quantities through the 
numerical domain. After the period of initial transients, the places where non-conservative 
effects can be found are almost exclusively restricted to the edge of the axial cut-out and 
the region roughly 45° from the polar axis within the ergosphere. 

We have verified that our new code is second-order accurate for smooth solutions and 



satisfactorily passes the tests described in lGammie et al.l (120031 ). A quantitative comparison 
of our code's performance to that of GRMHD will be left for future work. 
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3. Results 

Our initial condition is a torus of gas in hydrostatic equilibrium, entirely contained 
within the simulated volume; our goal is to present results characteristic of an accretion flow 
with a fixed aspect ratio in a long-term equilibrium. Before quoting results directly from the 
simulation data, we must therefore do two things: demonstrate that the fixed aspect ratio 
is achieved, and define more precisely the degree to which the simulation is in a statistical 
steady-state with respect to inflow. 

3.1. Scale-height regulation 

We set the parameters of our cooling function so that the ratio of the sound speed to 
the local orbital speed would produce a disk with a constant aspect ratio H/r = 0.13. In 
Figure [IJ we show how well the temperature was held to T* by comparing the time-averaged 
volume- weighted temperature in the bound accretion flow to the local value of T*. In the 
main disk body, this mean value was about (0.93-0.95)7^, but it rises sharply inside the 
ISCO. In other words, our cooling function succeeded in holding the disk temperature very 
close to (in fact, slightly below) the target temperature, but inside the ISCO, where the 
inflow time becomes comparable to or shorter than the cooling time, the temperature rises 
well above T*. 

How well our temperature-regulation led to a disk aspect ratio matching the goal value 
of 0.13 can be seen in Figure [2J The actual H/r was slightly above the goal (~ 0.14) through 
most of the simulation volume, but with a tendency to diminish inward inside r = 20M. At 
r = 10M, H/r ~ 0.12; by the time the flow reaches the ISCO, it is only ~ 0.07. Comparison 
with the curve showing how the scale-height changes as a result of including the relativistic 
correction to the vertical gravity (as discussed in § 12. 3p demonstrates that this thinning at 
small radius can be largely attributed to neglect of that effect. Thus, use of our cooling 
function achieved its principal goal: to place the scale-height of the disk under explicit 
control. 

Because our cooling function has a target temperature depending only on radius, at any 
particular radius the gas in the main body of the disk is nearly isothermal, and the density 
profile is therefore close to Gaussian (Fig. [3]). At higher altitudes above the midplane, the 
density falls slower than the Gaussian, presumably due to magnetic support. For this reason, 
the moment scale- height is slightly greater than the HWHM (Fig. [2]). 

We chose a value of H/r small enough that a key approximation of the NT theory 
could be approximately replicated in the simulation: the prompt radiation of dissipated heat. 
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Fig. 1. — Ratio of mean temperature (T — l)(e) to target temperature T*. The time-averaging 
interval was 7000-1 5000M. 
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Fig. 2. — Time-averaged density scale-height as a function of radius (solid curve), and time- 
averaged HWHM (dotted curve). The data were sampled every 20M from t = 7000M to 
15000M. Hydrostatic scale-height assuming the shell- and time-averaged temperature but 
employing the relativistic correction described in § I2.3I (dashed curve). 
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However, if the disk is to have a finite thickness, it cannot radiate all its heat. The parameters 
we chose for the cooling function yielded an accretion rate-weighted mean specific enthalpy 
that was well-described by h ~ 1 + 0.031(r/M)~ a8 . At large radius, where the Newtonian 
approximation applies, the ratio of h — 1 to the net binding energy is ~ 0.06(r/M) 2 , while 
at the ISCO this ratio is ~ 0.1. Thus, this toy-model does assure that the majority of the 
dissipated heat is radiated. 

3.2. Inflow equilibrium 

If the accretion flow were in a strict steady-state, the local (i.e., shell-integrated) mass 
accretion rate M(r) would be the same at all radii at all times and the mass interior to 
a given radius would likewise be constant. In these turbulent disks fed by a finite mass 
reservoir, the most we can hope for is that the time-average local accretion rate is nearly 
constant as a function of r through most of the accreting region, and the mass of the inner 
disk, after an initial period of growth, eventually levels off and fluctuates within some range. 
The degree to which we approach these goals is shown in Figures H] and [5j In the left-hand 
panel of Figure HI we see that the accretion rate (measured at the event horizon) varies by 
roughly a factor of five in an extremely irregular way. Nonetheless, as shown in the right- 
hand panel, the time-averaged M(r) is very nearly constant from the horizon to r ~ 14M for 
the latter 8000M of the simulation. The reason why we choose the interval 7000M-15000M 
for averaging is shown in Figure [51 As this figure demonstrates, it takes roughly the first 
7000M of the simulation for the mass of the inner disk to reach a rough plateau. Because 
the mass interior to a given radius fluctuates, we chose the starting point for time-averaged 
quantities to be the point at which essentially all the inner disk had reached at least 90% of 
its final mass, which is approximately t = 7000M. 

However, for the purposes of estimating the radiative efficiency, we require a tighter 
definition of inflow equilibrium. This is because we wish to contrast the computed radiation 
rate with the NT rate at an accuracy of a few percent or better. In the NT model, 23% 
of the total light is emitted between 12M and 25M, where our simulation shows significant 
departures from inflow equilibrium; a further 27% comes from outside 25M, where our 
simulation is not an accretion flow and we do not compute the luminosity at all. For these 
reasons, when we contrast the NT luminosity with that produced in the simulation, we 
adjust the local accretion rates to mimic inflow equilibrium and attach a carefully-chosen 
representation of large-radius emission where needed (see below for details). 
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Fig. 3. — Time- and azimuthally-averaged density (solid curves) at the ISCO (black) and 
r = 12 M (blue), each fit to a Gaussian (dashed curves). 
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Fig. 4. — (left) The accretion rate (in code units) at the event horizon as a function of time. 
A rate of 0.005 translates to accreting a fraction 0.14 of the initial mass in a time of 10000M. 
(right) The shell-integrated accretion rate as a function of radius, averaged from t = 7000M 
to 15000M, sampled every 1M. 
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3.3. Explicit radiative efficiency 

The fluid- frame emissivity Fff(r) found in the simulation and the NT prediction for 
this quantity are displayed in Figure [61 In the leftmost panel, we show how they com- 
pare when the NT emissivity uses the time-averaged accretion rate at the horizon over the 
same interval for which the simulation data were averaged, i.e., t = 7000M-15000M. For 
precise comparison of the two radiation models, we remove the effects of deviations from 
inflow equilibrium by altering the NT emissivity so that the value at any given radius corre- 
sponds to the time-averaged accretion rate at that radius, as determined by the simulation. 
Put another way, the fluid-frame surface brightness in a truly time-st eady NT model may 



be written as (3/47r)(GMM/r 3 ).R R (r) (notation as in IfCrolikl (jl999al )); we adjust this to 
(3/47r)(GMM(r)/r 3 ) J R iJ (r), with M(r) the time -averaged accretion rate at radius r in the 
simulation. By doing so, we compare the two radiation models in a way that factors out any 
contrasts due solely to fluctuations in the accretion rate. The adjusted version is shown in 
the right-hand panel of Figure El 

As this pair of figures shows, the two models coincide closely in the main disk body, but 
contrast sharply near and within the ISCO, which is at r m 2.3M for this spin. Because the 
NT model is founded on energy and angular momentum conservation in a time-steady disk, 
this coincidence is no surprise where the influence of the NT no-stress boundary condition is 
small. The principal departure between the two, a systematic offset in which the simulation 
curve lies ~ 10% below the NT curve, is due to the small fraction of the dissipated heat 
that the gas must retain to provide vertical pressure support in the disk. Near the ISCO, 
the accretion-weighted mean specific enthalpy is ~ 0.018 greater than unity. This thermal 
energy is 12% of the binding energy at the ISCO (0.155 per unit rest-mass). 

At small radius, however, the fluid-frame surface brightness of the simulation differs 
substantially from the NT model. At r = 3M, the simulation surface brightness is greater by 
40%; at the ISCO, although the NT model would predict no radiation, the surface brightness 
indicated by the simulation is roughly the same as at r = 3M; close to the horizon, the surface 
brightness rises to about twice the maximum predicted by the classical model. 

Another way to characterize the contrast between the simulation resul ts and the NT 



model is through the intermediary of another analytic model. In the model of lAgol fc Krolik 



(120001 ) . it is supposed that a finite stress is exerted at the ISCO, but all other assumptions 
follow those of NT. This model (which we will abbreviate as AK) is parameterized by the 
additional efficiency Ae due to the non-zero stress at the ISCO; in the curves shown in the 
two panels of FigureEl Ae = 0.01, a value chosen as an approximate best fit between the AK 
model and the simulation data. In the region immediately outside the ISCO, where the AK 
model is defined, it does a reasonable job of reproducing the simulation results, particularly 
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Fig. 5. — The mass contained within four sample radii: 14M (solid curve), 12M (dotted 
curve), 10M (dashed curve), and 8M (dash-dot curve), all as functions of time. The thin 
solid lines mark 90% of the final mass for each of these radii. A mass of 10 in code units is 
2.8% of the initial torus mass. 




Fig. 6. — Radiated flux per unit area in the fluid frame as a function of radius: time-averaged 
simulation data (solid curve); as predicted by the NT model (dotted curve); as predicted 
by the AK model with Ae = 0.01 (dashed curve), (left) Using the time-averaged data from 
7000M to 15000M. (right) Adjusting the NT and AK emissivities as described in the text. 
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when allowance is made for the retained heat. 

Only some of this radiation reaches infinity, and any that does arrives with a significant 
Doppler shift, most often toward the red. Using the techniques described in § 12.31 and 
Appendix A, we computed the luminosity received at infinity per unit radial coordinate 
dL/dr, which is shown in Figure [7J Like the emissivity in the fluid frame, dL/dr for the 
simulation data in the main disk body closely tracks the NT prediction. The only difference 
between the two is that the simulation data version lies slightly (~ 10%) below the NT 
curve: this offset is simply another reflection of the offset already seen in the fluid-frame 
emissivity due to the non-zero heat content of a physical disk. At small radii, the shelf in 
the fluid-frame emissivity is transformed into an inward extension of significant luminosity 
that extends from r ~ AM to r ~ 2M. Although the fluid-frame emissivity extends farther 
inward, its efficiency in creating luminosity at infinity is cut off by a combination of increasing 
redshift and probability of photon capture by the black hole. 

At larger radii, departures from inflow equilibrium become significant. To compute 
accurately J dr dL/dr, the data of our simulation must be both adjusted so as to correspond 
to true inflow equilibrium and supplemented by an extension to larger radius to account for 
the substantial radiation from radii larger than 25M. Because the time-averaged fluid-frame 
emissivity in the simulation tracks the NT model so closely for 5M < r < 12M, we define 
the simulation luminosity as its dL/dr integrated from the horizon to r = 12M plus the NT 
luminosity at the mean accretion rate from r = 12M outward. 

Given that definition, we find that the efficiency with which this simulation generated 
light reaching infinity, averaged from 7000-15000M, was 0.151. This number is 6% greater 
than the NT figure, which is 0.143 after allowing for photon capture. 

3.4. Extrapolating to the complete radiation limit 

As discussed in the Introduction, our principal goal in this initial simulation was to 
achieve as close a test as possible of the effect of the ISCO stress boundary condition on 
the radiative efficiency. We must now evaluate the degree to which our only approximate 
replication of the other NT assumptions affected this test. Time- and azimuthal- averaging 
should provide a good approximation to a stationary state and axisymmetry; incomplete 
radiation of the dissipated heat is our principal concern here. 

We have already seen that our radiation rate closely tracks the NT radiation rate in 
the main disk body, but is about 10% lower. Thus, to extrapolate to complete radiation, 
the emission from this portion of the flow should be increased by this amount. Near the 
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ISCO, where the effects of the stress boundary condition become important, we cannot use 
this comparison method to estimate the magnitude of the retained heat. Instead, we observe 
first that at the ISCO the mean Thomson optical depth through the disk in our simulation 
was ~ 500m, where m is the accretion rate in Eddington units. The corresponding diffusion 
time is ~ 0.7m orbits. At the same place, the inflow rate is ~ 0.6fi = 1.27r/P rb- Thus, the 
photon diffusion time near the ISCO in a real disk should be shorter than the inflow time — 
and shorter than our toy-model cooling time — for all accretion rates below Eddington. A 
second standard of comparison may be derived from the magnitude of the retained heat. We 
found earlier that the accretion-weighted mean specific enthalpy is ~ 1 + 0.018 at r ~ 2M. 
That the retained heat is ~ 10% of the binding energy there is consistent with the fact that 
~ 10% of the heat dissipated in the main disk body is left unradiated. Combining these two 
arguments, we might expect that in the limit of truly complete radiation of dissipated heat, 
the efficiency could have been greater by as much as 0.02, rising perhaps to ~ 0.17, 20% 
above the classical number as adjusted for photon capture. 

Additional heat is created in the plunging region (the mean accreted specific enthalpy 
rises from 1.02 at the ISCO to ~ 1.03 at the horizon), but, as we have already seen, the 
fraction of photons escaping from regions so close to the horizon to infinity is relatively small, 
so only a small part of the additional 0.01 in rest-mass equivalent is likely to reach distant 
observers. 

We might also ask what effect truly radiating all the heat would have on electromagnetic 
energy fluxes. To approach this question we begin by considering it from the point of view of 
the classical (NT) theory of accretion, where much attention is paid to the r-0 component of 
the stress tensor Tjf, but little is said about other components except for the assumpt i on that 



the stress tensor is orthogonal to the four- velocity, u^Tjf = 0. As iBeckwith et al.l (j2008al ) 
pointed out, this assumption is consistent with the sort of stress NT had in mind, i.e., or- 
dinary viscosity, but not necessarily with other physical stress mechanisms. In particular, 
it is inconsistent with MRI-driven MHD turbulence: the electromagnetic stress tensor con- 
tains a term ||6|| 2 w' 1 u l/ , which is manifestly not orthogonal to the four-velocity; in addition, 
the turbulence entails another (generally rather smaller) contribution to the stress tensor 
(ph+ WbW^Su^Siiv, where Su^ is the fluctuating part of the four-velocity. Described in more 
qualitative terms, the classical theory accounts for the energy flow due to the work done by 
the stress, but not the energy flow due to the advection, by the mean flow, of an energy 
density associated with the stress mechanism. 

As numerous numerical studies of the MRI-driven turbulence have shown, the fluid- 
frame ratio a mag = 2(b r b < t > )/(\ \b\ I 2 ) ~ 0.2-0. 3 in th e disk body, rising by factors of a few in 
the plunging region (e.g. Hawley fc Krolik ( 20021 )). At the order of magnitude level, the 



- 21 - 



ratio of the advected magnetic energy flux to the magnetic work is ~ u r / (a mSug ru' l) ) , which 
is very small in the disk body, but rises sharply near the ISCO and in the plunging region. 
In this simulation, we find that the time-averaged advected magnetic energy flux per unit 
rest-mass is 0.03 at the ISCO, a significant contribution to the energy budget. 

To complete our extrapolation to complete radiation therefore means that we need to 
determine how u r near the ISCO might change when that limit is taken at fixed accretion 
rate. Fixing the accretion rate means that the vertically-integrated stress does not change, 
and we do not need to estimate how the stress would change as a function of the disk's 
thermal state. If u r near the ISCO depends primarily on the shape of the potential, the 
advected magnetic energy flux per accreted rest-mass would remain roughly the same. On 
the other hand, if u r in this region depends on the gas thermal content in the sense that it 
increases with increasing temperature, more complete radiation would also lead to a smaller 
rate of magnetic energy advection, and therefore to a larger net outward Poynting flux and 
a larger amount of energy available for dissipation. Which of these possibilities lies closer to 
the truth (and under which circumstances) remains to be determined. 



4. Summary and Implications 

Global disk simulations have for many years focused on dynamical effects, i.e., angular 
momentum transport leading to inflow. To link them to observations, however, requires 
including considerations of thermodynamics, for the energy to radiate photons is drawn 
from the thermal energy of the gas (whether or not the particle distribution functions are in 
fact near those of thermal equilibrium). By combining an energy conserving algorithm with 
an explicit cooling function in a new simulation code, HARM3D, we are able to begin the first 
steps toward drawing that connection. 

In this first application of our new technique, we have found that a disk with H/r ~ 0.1 
accreting onto a black hole with spin parameter a/M = 0.9 carries thermal and magnetic 
energy past the ISCO at a rate ~ 0.05 per unit rest-mass, while producing radiation that 
reaches infinity at a rate ~ 0.15 per unit rest-mass. These numbers contrast with those of 
the classical NT model, in which the flow carries no thermal or magnetic energy, and for 
a/M = 0.9 radiates ~ 0.14 per unit rest-mass to infinity. Determining the observed luminous 
efficiency of a more realistic accretion disk, as opposed to the ideal NT model, depends on 
the careful assessment of several potentially offsetting effects. First, additional thermal, 
magnetic, and radiated energy can be drawn from the orbital energy by magnetic stresses 
that can persist through the location of the ISCO and all the way down to the horizon. 
However, only a fraction of that energy need be radiated, with much of the remainder 
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retained as heat and magnetic field captured by the hole. Next, even if there is enhanced 
photon production near and inside the ISCO, for this particular spin, the combination of 
comparatively high capture probability and gravitational redshift means that little radiation 
from inside the ISCO reaches infinity For lower spin holes the ISCO is further from the 
horizon and the plunging regi on can be more effectively represented in the luminosity at 
infinity faeckwith et al.ll2008ah . 



These results have implications for the spectral shape of the emitted radiation. Generi- 
cally, the effect of the continuing stresses is to move the radius of peak emission inward and 
raise the fluid-frame effective temperature at that location. For example, in this instance the 
maximum in dL/dr (after allowing for photon capture and all Doppler shifts) moves from 
the NT prediction of r ~ 4.3M to r ~ 3.5M. Similarly, the fluid-frame flux at the peak of 
dL/dr is about 30% greater (7% higher effective temperature) in the simul ation data than in 
the N T mo del. In terms of the rad iation edge terminology introduced by iKrolik fc Hawley 
(120021 ) and iBeckwith et al.l (j2008al ). we find that 95% of the radiation reaching infinity is 
produced outside r = 2.75M, in contrast to 3.6M in the NT model. 



In a previous study, IBeckwith et al.l (j2008al ) used the stress distributions observed in an 
ensemble of disk simulations to estimate the dissipation that might be associated with those 
stresses, and from this the accretion efficiency and maximum temperature in the spectrum 
reaching infinity. After accounting for photon capture and Doppler-shifting effects, they 
found that, depending on the particular simulation examined and the topology of the initial 
magnetic field, the luminosity reaching infinity could be anywhere from 20%-100% greater 
than NT when a/M = 0.9. The low end of this range was produced by an accretion flow 
whose initial field was entirely toroidal, the high end by an accretion flow whose initial field 
was a large dipolar loop, as in the present simulation. Thus, there is a sizable gap between 
their estimate of the radiative efficiency and ours. 

Applying Beckwith et al.'s expression to our data leads to a prediction for the dissipation 
rate very similar to theirs^. The fact that our radiation rate is considerably less than this 
prediction suggests that the simple ansatz used by Beckwith et al. to directly compute 
dissipation from stress and equate dissipation with radiation is a simplification that likely 
overestimates the net emission. For example, because our cooling function's radiation rate is 
at most comparable to the inflow rate near and inside the ISCO, not all the heat dissipated 
in that region can be radiated. However, even if all the heat generated in this simulation 
were radiated, the increase in efficiency relative to NT would be only ~ 20%. In addition, 



7 In fact, the accretion rate histories of the two simulations are remarkably similar, suggesting that the 
underlying physics imposes a long-term order despite the significant difference in computational algorithms. 
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not all the work done by the stress necessarily goes into a form that is dissipated. Kinetic 
and magnetic energies can be advected with the accretion flow into the black hole, producing 
no effective increase in overall efficiency. This point is closely related to the issue of advected 
energy discussed in § 13.41 

Framed in the context of predictions for real accretion flows in Nature, these questions 
emphasize the importance of realistic dissipation and radiation physics for obtaining more 
accurate accounts of radiation associated with accretion. In the vicinity of the ISCO, where 
the energy available for release is largest, one cannot say with confidence that in general the 
dissipation and cooling times are shorter than the inflow time. Moreover, both processes are 
likely to depend on the detailed circumstances pertaining to any particular accreting black 
hole, so that there may not be a single efficiency number applicable to all black holes of a 
given spin. 

In sum, we have shown that by use of a toy-model optically-thin cooling function, it 
is possible both to control the thickness of the accretion flow and to tally (approximately) 
the rate at which radiation can be produced by dissipation in the flow. At relatively large 
radii, where the inflow time is long compared to the cooling time, our ansatz of substituting 
gridscale dissipation for genuine microphysics and radiating the heat so generated at an 
arbitrarily chosen rate is capable of capturing the global energetics of accretion reasonably 
well. However, at smaller radii (particularly near and inside the ISCO), where the inflow 
time can be comparable to the cooling time, use of realistic dissipation and radiation rates 
can be more important. 

Having demonstrated the technical feasibility of this approach, we will next employ 
it to explore more fully how accretion onto black holes depends on disk thickness and on 
black hole rotation rate. In this context, we point out that although there is a standard 
notation for describing black hole rotation (the spin parameter a/M), there are several 
extant definitions of the scale- height, differing from one another by factors of order unity. 
We use the vertical density moment; standardization of this definition would be of benefit 
so that different calculations can be compared quantitatively without confusion. 

Lastly, we remark that in this paper we have set aside the fact that photons are not the 
only form in which energy can be sent to infinity from the vicinity of black holes. Accreting 
black holes are also capable of driving mass motions, often relativistic, that can carry sig- 
nificant power in Poynting flux. Simulational work exploring the a s sociated luminosity ha s 



already begun (IMcKinnev fc Gammiell2004l ; lHawlev fc Krolikl 120061 : iBeckwith et al.ll2008bh . 



In future work, we will use the new simulation code introduced here to relate the energetics 
of those outflows more closely to the accretion energy budget. 
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A. The Radiative Transfer Calculation 

Our me t hod f or calculating the radiative transfer closely follows the one described by 



Noble et al.l (120071 ). We have made many changes to the code, including the ability to use 
3D simulation data, different emission models (such as the one explained here), and many 
optimizations that have made it significantly faster. 

The algorithm integrates geodesies from the observer's camera through the source 
domain — our simulation data. These geodesies point toward the camera and the future. 
A geodesic represents a path along which a bundle of photons travel. The Lagrangian form 
of the geodesic equations is used: 

dx^ 8N 

^ = N" > it^nNuW , (Al) 

where x^ is the world-line of the photon bundle and is the geodesic's tangent vector 
parameterized by the affine parameter A. 



This preprint was prepared with the AAS IATgX macros v5.2. 



Since there is no absorption or scattering, the radiative transfer equation takes the form 



dl 



(A2) 



where X = I u /v z is the Lorentz invariant intensity, l v is the specific intensity, J = jy/v 2 
is the invariant emissivity, j u is the emissivity, and v is the local frequency of the photon. 
For the purposes of calculating the bolometric luminosity, we consider only line emission. 
We can assume either constant emission frequency (e.g. Fe Ka fluorescence) where we must 
integrate over all frequencies at the camera, or constant observer frequency where we assume 
the emission is contrived to emit at a frequency which — when redshifted to the camera's 
frame — is equal to the frequency of observation. It is easy to show that both methods give 
the same bolometric luminosity. We therefore choose the latter method as it requires less 
computational effort. 

We assume that the radiation is emitted isotropically, so j v oc £/ (4tt) but we must 
also take into account the constraint that the fluid's emission frequency is the blueshifted 
frequency at the observer: 



where G(X), the redshift factor, is the ratio of the photon's energy measured by the camera 
to the photon's energy measured by the fluid: 



4ttv 2 



5{v-v /G{\)) 



(A3) 





(A4) 



Here, is velocity of the camera which is assumed to be static in flatspace; this is a good 
approximation as we place the camera 10 6 M away from the black hole. 
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Fig. 7. — Luminosity received at infinity per unit radial coordinate: time-averaged simulation 
data (solid curve); as predicted by the NT model (dotted curve), (left) Using the time- 
averaged data from 7000M to 15000M. (right) Adjusting the NT emissivity as described in 
the text. 



